
addpath './Model/'
addpath './Visualization/'
addpath './Util/'
addpath './BoundaryIntegral/'
addpath './KernelGradient/'

global shiftAlpha shearModulus poissonRatio mesh;

shearModulus = 1;

poissonRatio = 0;

shiftAlpha = 1e-3;


mesh = readMesh('cube3x3.obj');

rc = mesh.v;

rx = rc(:, 1);
ry = rc(:, 2);
rz = rc(:, 3);

echo = 'mesh loaded'
% 
A = assemblyMatrix(mesh);

echo = 'matrix assemblied'
% 
[b] = assemblyBoundaryVector(mesh); 

echo = 'boundary assemblied'
% 
theta = A \ b;
% 
% theta = bicgstab(@matrixVectorProduct, b, 1e-4, 100);
% 
% % theta = gmres(@matrixVectorProduct, b, 10, 1e-4, 100);
% 
echo = 'matrix solved'
% 
u = evaluateInteriorPoint(rx, ry, rz, mesh, theta);

echo = 'interior points evaluated'
% 
plotDeformation(mesh, rx, ry, rz, u);

